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A.  ACCOMPLISHMENTS  TO  DATE 

The  following  is  a  brief  summary  of  accomplishments  to  date  under  Grant 
AFOSR-84-0131. 

1.  The  Reduced  Basis  Method  for  Algebraic  Systems 

'  \ 

1.1  Subspace -Projector  Pairings .  Implementation  of  the  reduced  basis  method 
requires  the  choice  of  a  subspace  and  a  projector  onto  that  subspace.  For  an 
arbitrarily  chosen  subspace -projector  pair,  existence  of  the  true  solution 
curve  is  not  sufficient  to  guarantee  the  existence  of  the  corresponding 
reduced  basis  solution  curve.  However,  when  the  former  curve  exists,  it  has 
been  shown  in  fAl]  that  there  are  infinitely  many  subspace -projector  pairings, 
each  utilizing  an  arbitrarily  selected  subspace ,  under  which  the  reduced  basis 
solution  curve  exists.  Moreover,  the  resulting  error  estimates  are  of  the 
same  nature  as  those  that  apply  in  the  more  familiar  case  when  a  subspace  is 
paired  with  its  orthogonal  projector.  __ 

1.2  Reduced  Basis  Additive  Correction  Methods.  Additive  Correction  Methods 
have  been  considered  by  a  number  of  authors  as  a  means  of  accelerating  slowly 
converging  iterative  processes  (see,  for  example,  [A2]-[A4]).  Furthermore,  it 
has  been  recognized  that  additive  correction  is  central  to  the  basic  idea  of 
multigrid  methods  [A4-A5].  Although  the  reduced  basis  method  in  its  original 
form  appears  to  have  little  in  common  with  additive  correction,  a  class  of 
such  methods  has  been  developed  using  the  reduced  basis  concept  [1].  Further¬ 
more,  it  has  been  shown  that  in  their  "two-grid"  form,  certain  multigrid 


methods  are  special  cases  of  this  class.  The  reduced  basis  point  of  view  pro¬ 
vides  insight  into  the  error  reduction  capability  of  such  multigrid  methods 
and  suggests  additive  correction  variants  that  may  be  more  effective  than 
those  commonly  employed  in  multigrid  methods. 

2.  The  Reduced  Basis  Method  for  Systems  of  Differential  Equations 

Error  Estimates  for  the  reduced  basis  method  solution  of  differential  and 
differential-algebraic  equation  systems  are  contained  in  the  thesis  [2]. 
These  estimates  are  local  in  the  following  two  senses.  First,  they  apply  on  a 
nontrivial,  t  it  possibly  very  small  interval.  Second,  they  require  some  point 
of  the  reduced  curve  to  lie  on  the  true  solution  curve.  The  recent  research 
reported  in  [18]  has  removed  the  interval  length  restriction  in  the  differen¬ 
tial  equation  case  and  extended  the  error  analysis  to  global  versions  of  the 
methods  in  [2],  thus  effectively  eliminating  the  second  local  aspect  of  that 
work.  Furthermore,  this  work  also  incorporates  the  effects  of  the  errors 
resulting  from  the  numerical  integration  of  the  reduced  basis  ODE  systems. 

3.  Two -Fluid.  Two -Phase  Flow 

Additional  theoretical  results  on  the  nature  of  the  void  fraction  have 
been  incorporated  into  [A6]  resulting  in  the  revision  [3]. 

4.  Binary  Gas  Mixture  Flow  Through  Combustors 

In  an  attempt  to  reduce  the  development  cycle  costs  associated  with 
design  of  gas  turbine  engine  combustors,  mathematical  combustor  models  are 


being  employed  to  provide  information  about  performance  trends  and  to  predict 
velocity,  pressure  and  thermodynamic  property  profiles  in  simulated  practical 
combustion  environments.  It  has  been  demonstrated  that  the  dual  variable 
method  can  be  applied  to  the  predictive  model  of  the  fluid  dynamics  associated 
with  an  axially  symmetric  centerbody  combustor  being  studied  at  WPAFB.  This 
work  appeared  in  [4] . 

5 .  Error  Estimation  and  Singularities 

The  central  theme  of  this  project  concerns  the  discretization  error  for 
general  nonlinear,  parameter -dependent  equations  of  the  form  F(z,A)  -  0  where 
F  is  a  nonlinear  mapping,  z  is  the  state  variable  representing  the  solu¬ 
tion,  and  A  is  a  vector  of  parameters  that  characterize  intrinsic  properties 
of  the  system  or  extrinsic  quantities  influencing  its  behavior.  In  the  case 
of  fluid  problems,  the  operator  F  may  be  generated,  for  example,  by  the 
time -independent  Navier  Stokes  equations  together  with  the  necessary  boundary 
conditions . 

In  general,  the  solution  set  of  such  parametrized  equations  constitutes  a 
differentiable  manifold  of  dimension  equal  to  that  of  the  parameter  space. 
While  this  fact  is,  of  course,  well  known,  we  appear  to  have  been  the  first 
who  have  been  using  this  fact  as  the  basis  for  a  consistent  study  of  the 
numerical  problems  for  these  equations.  Our  results  have  begun  to  show 
clearly  the  value  and  power  of  this  geometric  approach. 

For  numerical  considerations,  an  important  aspect  of  this  theory  concerns 
the  estimation  of  the  "distance"  between  the  manifolds  defined  by  a  given  dif¬ 
ferential  equation  and  by  its  discretization,  respectively.  This  has  been  the 


topic  of  a  series  of  papers  by  J .  P.  Fink  and  W.  C.  Rheinboldt  with  partial 
support  under  this  grant.  See  the  earlier  articles  [A7],  [A8]  and  then  [6] 
and  [16].  In  particular,  in  the  last-named  paper  [16]  we  have  been  the  first 
to  analyze  the  case  of  multi -dimensional  manifolds  which  is  of  increasing 
importance  in  applications. 

An  essential  aspect  of  these  studies  concerns  the  question  as  to  the 
exact  definition  of  the  error  between  a  solution  manifold  and  its  approxima¬ 
tion.  Obviously,  the  error  depends  on  which  points  are  to  be  compared.  In 
the  cited  papers  it  was  shown  that  this  correspondence  between  the  points  on 
the  two  manifolds  has  to  be  defined  by  appropriate  local  coordinate  systems. 
In  other  words,  the  resulting  error  is  controlled  by  the  choice  of  the  local 
coordinate  system,  and,  since  the  error  measure  must  be  physically  meaningful, 
not  all  local  coordinate  systems  are  equally  desirable.  This  question  becomes 
particularly  acute  in  the  vicinity  of  singular  points  where  the  behavior  of 
the  solutions  may  be  subject  to  change. 

This  connection  between  the  choice  of  the  local  coordinate  systems  and 
the  singularity  behavior  of  a  point  has  led  us  to  a  closer  study  of  admissible 
coordinate  systems  at  foldpoints.  But  the  results  in  the  papers  [5]  and  [7] 
also  suggested  that  the  open  questions  about  the  proper  choice  of  the  coordi¬ 
nate  system  for  physically  meaningful  error  measures  requires  a  much  closer 
look  at  the  characterization  of  foldpoints  on  manifolds  and  at  the  methods  for 
their  computations.  This  is  the  topic  of  the  next  subsection. 


Detection  and  Computation  of  Singularities 


In  addition  to  its  application  to  error  estimation,  singularity  theory  is 
also  being  used  to  develop  methods  for  the  detection  and  computation  of 


singular  points  on  the  solution  manifold  of  a  nonlinear  parameter -dependent 
equation  F(z,A)  -  0. 

As  noted  in  the  previous  subsection,  the  solution  set  of  a  parametrized 
equation  F(z,A)  -  0  represents,  in  general,  a  differentiable  manifold  in  the 
combined  space  of  the  state  variable  and  the  parameter  vector.  This  requires 
a  regularity  assumption  which  is  not  very  restrictive  in  applications,  but 
which  --  from  the  viewpoint  of  singularity  theory  --  implies  the  use  of  suit¬ 
able  "unfoldings" .  In  most  practical  problems  a  host  of  very  natural  unfold¬ 
ings  suggest  themselves.  The  particular  choice  of  unfolding  parameters 
affects  the  location  and  type  of  the  resulting  fold-points  on  the  manifold  and 
with  it  also  the  questions  raised  in  the  previous  subsection. 

In  the  past  years,  our  work  in  this  area  concentrated  first  on  the  compu¬ 
tation  of  foldpoints  by  means  of  one  of  the  many  possible  forms  of  augmented 
equations.  One  approach  in  this  direction  was  the  use  of  the  tangent  map  of 
differential  geometry  which  was  exploited  in  the  already  mentioned  papers  [5] 
and  [7]  and  led  there  to  a  geometrically  Instructive  and  coordinate  free 
treatment  of  fold  points,  in  general. 

But  these  results  also  pointed,  once  again,  at  the  need  for  a  more 
detailed  study  of  the  geometrical  aspects  of  the  singular  points.  In  [21]  it 
was  shown  that  it  is  possible  to  reformulate  in  our  setting  some  of  the  basic 
results  of  bifurcation  theory  related  to  computational  aspects.  The  availa¬ 
bility  of  solution  manifolds  opens  up  the  various  tools  of  differential 
geometry  and  provides  a  very  natural  framework  for  the  desired  reformulation 
of  the  theory.  In  particular,  the  bifurcation  sets  appear  as  cuts  of  the 
solution  manifold.  This,  in  turn,  corresponds  again  to  the  consideration  of 
augmented  equations  and  hence,  our  new  theory,  does  provide  indeed  a  new 
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approach  to  the  vexing  problem  of  the  properties  of  the  numerous  augumentation 
considered  in  the  literature.  In  [21]  the  results  were  already  used  to 
analyze  a  particularly  promising  augmentation  for  the  computation  of  an  entire 
class  of  fold  points.  This  work  is  now  continuing.  Some  further  results  on 
this  topic  were  presented  in  [19]. 

The  general  computational  problem  in  this  area  is,  of  course,  the  deter¬ 
mination  of  the  principal  characteristic  features  of  the  solution  manifold  of 
the  given  parametrized  equation.  This  includes,  in  particular,  the  computa¬ 
tion  of  the  foldpoints,  but  covers  also  other  features.  The  principal  methods 
used  today  for  analyzing  such  manifolds  are  the  continuation  methods.  But 
before  any  such  method  can  be  applied  in  the  case  of  a  multi -dimensional  mani¬ 
fold,  the  parameter  dimension  has  to  be  reduced  to  one.  Geometrically,  such  a 
reduction  is  equivalent  with  a  restriction  to  some  path  on  the  manifold  of  the 
original  equation.  A  continuation  method  then  computes  a  sequence  of  points 
along  such  a  path.  Clearly,  in  general,  it  is  not  easy  to  develop  a  good  pic¬ 
ture  of  a  p-dimensional  manifold  from  information  along  one -dimensional  paths. 
Thus,  it  is  not  surprising  that  there  is  growing  interest  in  computational 
methods  which  generate  multi-dimensional  grids  of  solution  points  covering  an 
entire  portion  of  the  manifold.  In  [22]  a  new  method  for  this  purpose  was 
presented.  It  allows  for  the  computation  of  the  vertices  of  a  simplicial  tri- 
angulation  of  a  p-dimensional  solution  manifold  of  a  parametrized  equation. 
An  essential  part  of  the  method  is  a  constructive  algorithm  for  computing  mov¬ 
ing  frames;  that  is,  of  orthonormal  bases  of  the  tangent  spaces  that  vary 
smoothly  with  their  points  of  contact.  The  triangulation  algorithm  uses  these 
bases,  together  with  a  chord  form  of  the  Gauss-Newton  process  as  corrector,  to 
compute  the  desired  vertices.  The  Jacobian  matrix  of  the  mapping  is  not 
required  at  all  the  vertices  but  only  at  the  centers  of  certain  local 
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"triangulation  patches".  Numerical  experiments  have  already  shown  that  the 
method  is  very  efficient,  even  around  singularities.  This  opens  up  new  possi¬ 
bilities  for  determining  the  form  and  special  features  of  such  solution  mani¬ 
folds  . 

7.  Finite  Element  Formulation  of  the  Streamfunction-Vorticity  Equations 

The  Navier- Stokes  equations  can  be  written  in  primitive  variable  formula¬ 
tion,  in  terms  of  the  streamfunction  as  a  fourth  order  problem  or  as  two 
second  order  equations  in  the  streamfunction-vorticity  formulation.  In  the 
linear  case  the  fourth  order  problem  for  the  streamfunction  is  the  well-known 
biharmonic  equation.  Although  the  primitive  variable  formulation  has  received 
the  most  attention,  the  streamfunction-vorticity  formulation  is  also  of  con¬ 
siderable  interest  in  two  dimensional  domains.  That  is  partly  due  to  the  fact 
that  only  two,  as  opposed  to  three,  fields  are  to  be  computed;  but  it  is 
mainly  due  to  the  fact  that  the  incompressibility  constraint  Is  avoided 
through  the  introduction  of  the  streamfunction. 

Several  theoretical  and  practical  issues  arising  in  the  finite  element 
approximation  of  the  streamfunction-vorticity  equations  have  been  studied. 
Initially  error  estimates  for  the  linear  problem  were  investigated.  Since  the 
velocity  is  expressed  in  terms  of  the  derivatives  of  the  streamfunction,  it  is 
of  practical  concern  to  ascertain  if  these  derivatives  are  optimally  approxi¬ 
mated  for  choices  of  elements.  Previous  analyses  concerning  this  problem  were 
improved  upon  and  the  optimality  of  the  error  verified  in  [A10] . 

Other  issues  arising  in  the  finite  element  approximation  of  the  stream- 
function  and  vorticity  include  computations  in  multiply  connected  domains,  the 
use  of  low  order  elements,  the  incorporation  of  a  variety  of  boundary 
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conditions  into  the  weak  formulation,  estimates  for  the  errors  in  the  finite 
approximations  for  the  nonlinear  problem  and  the  recovery  of  the  primitive 
variables.  A  preliminary  report  on  computations  in  multiply  connected  domains 
using  low  continuity  finite  element  spaces  was  presented  in  [11].  A 
comprehensive  report  dealing  with  all  of  the  theoretical  and  practical  issues 
mentioned  above  as  well  as  presenting  numerical  examples  is  given  in  [12]. 


8 .  A  Finite  Element  Analysis  of  MHD  Flow 

The  equations  governing  the  steady  flow  of  incompressible  electrically 
conducting  fluids  in  the  presence  of  a  given  magnetic  field  can  be  expressed 
as 


-  Au  +  |  (uV  u)  +  -  (BxV^)  -  (uxBxB)  -  0 

n 

-L4>  +  div(uxB)  -  0 

div  u  -  0 

where  u  is  the  velocity,  p  the  pressure,  ^  the  electric  potential,  B  the  mag¬ 
netic  field  and  N,M  given  dimensionless  parameters.  By  rewriting  certain 
terms  using  vector  identities  and  using  appropriate  spaces,  one  can  obtain  a 
weak  formulation  for  this  problem  that  is  similar  to  one  for  the  Navier-Stokes 
equations  written  in  terms  of  primitive  variables  (see  [All]).  The  purpose  of 
using  such  a  weak  formulation  is  to  take  advantage  of  the  results  already 

proved  in  [All],  Specifically,  the  weak  formulation  is  to  find  u  -  (u ,4>)  e  W, 
2 

p  t  L  such  that 
r  o 

a(u,v)  +  a.(u,u,v)  +  b(v,p)  -  (f,v)  for  all  v  t  W 


-  9  - 


b(u,x)  -  0  for  all  *  e  L 


where 


i(u,v)  -  -^  [  Vu:Vv  +  |  { V^-(uxB)  }{V\£>-(yxB) ) 

n  *  •* 

i1(u,u,v)  -  ^  J  {uV  uy  -  uV  vu} 


b(v,x) 


--J 


X  div  v 


112  2 
and  W  -  H  X  H  ,  L  -  {<f>  e  L  : 
—  — o  O  o 


■t- 


0). 


The  continuity  and  stability  conditions  necessary  to  guarantee  existence 
and  uniqueness  of  the  solution  of  the  weak  problem  have  been  proved.  In  addi¬ 
tion,  an  error  estimate  for  the  finite  element  approximation  of  the  weak  prob¬ 

lem  has  been  obtained.  These  results,  as  well  as  a  discussion  of  an  iterative 
method  for  solving  the  discrete  problem  will  be  presented  in  [13], 

9 .  Dual  Variable  Transformations 

The  dual  variable  method,  originally  developed  in  the  context  of  finite 
difference  discretizations  of  transient  incompressible  Navier-Stokes  equations 
[A9],  is  a  technique  to  considerably  reduce  the  size  of  the  linear  system  to 
be  solved  at  each  time  step.  The  method  involves 

(1)  the  determination  of  the  rank  of  the  discrete  divergence  operator, 
A, 

(2)  the  determination  of  a  basis  for  the  null  space  of  A,  and 

(3)  the  calculation  of  a  particular  solution  of  the  discrete  continuity 

equation. 

In  [8]  a  finite  element  implementation  of  the  dual  variable  method  is 
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presented  using  quadrilateral  piecewise  bilinear  velocity/constant  pressure 
elements.  Algorithms  for  the  determination  of  a  basis  for  the  null  space  of 
the  discrete  divergence  operator  and  a  particular  solution  are  presented. 

In  [9]  a  finite  difference  discretization  of  the  Navier- Stokes  equations 
describing  a  compressible  flow  problem  is  viewed  as  a  system  defining  flows  on 
an  associated  network.  This  observation  then  provides  a  means  of  applying  the 
dual  variable  method  to  economize  on  their  numerical  solution. 

The  nature  of  the  aerodynamics  in  and  around  such  structures  as  cavities 
and  deflectors  or  spoilers  on  various  aircraft  configurations  was  investigated 
using  the  dual  variable  method  [10]. 

A  summary  of  the  dual  variable  method  is  given  in  [14]. 

Iterative  methods  are  under  investigation  for  the  solution  of  the  linear¬ 
ized  finite  difference  discretizations  involved  in  the  dual  variable  method. 
The  generic  form  of  dual  variable  system  suggests  a  splitting  in  which  a 
Stieltjes  matrix  is  to  be  solved  at  each  step.  The  method  has  been  imple¬ 
mented  for  two  dimensional  domains  and  convergence  properties  are  being  inves¬ 
tigated  as  part  of  a  Ph.D.  thesis. 

10.  Fluid  Flow  on  Curved  Domains 

A  finite  difference  scheme  was  derived  for  two-dimensional,  transient, 

incompressible  Navier - Stokes  problems  in  which  the  flow  domain  Q  is  a  bounded 

2 

s i mpl v - connec ted  region  for  which  there  exists  a  C  invertible  mapping  r  onto 

t  he  uni t  square : 


r 


n  -  s  -  [0 , l ]x [ o , l ] 


y  «  l»  v*1  »T»-\,v;*-:r.-y;y^7’jy.-v.  y.  -  7 r*.  -it  r r. -7  7* 'FT  7?  -.»  ~7  '  ■», v 

-  11  - 


The  transformed  Navier-Stokes  problem  is 


div  l  •  4 

,ri 


-  0 


<f> ,  • 


/ 


,,t  +  Ijf  Vi/rj  “  -Pr^j.x,  +  -  </Uvu  .  )_  +  f,  i 


i/i 


|J|  VPJk  k.r^rj  'r  *i 


-  1.2 


subject  to  ^hitial  condition 

/ 


^  v(r ,0)  -  a(x(r) )  ,  r  t  S 

and  boundary  condition 


v(r,t)  -  0  ,  r  e  8  S  and  t  >  0  , 
where,  v(r,t)  is  velocity,  p  is  pressure,  the  Jacobian  matrix 


J  "  [xi  r 

’  J 


-  |J|J'  V  V 


±  -  |J|J  \  • 


The  finite  difference  discretization  of  the  above  equations  is  proven  to 
be  unconditionally  stable  and  convergent.  They  also  reduce  to  the  well  known 
Krzhivitski  and  Ladyzhenskaya  scheme  [A12]  for  rectangular  domains,  e.g.,  r 
the  identity.  This  work  was  the  subject  of  a  Ph.D.  Dissertation  [A13]  and  a 
recent  technical  report  [17]. 


11.  Equilibrium  Manifolds  in  Computational  Mechanics 


Equilibrium  problems  arise  naturally  in  continuum  and  fluid  mechanics  as 
a  set  of  nonlinear  equations  of  the  form 


F(z , \)  -  0 


(1) 


$ 

a 


£ 


£ 


a 


\  S  SV  VS 


where  z  is  from  the  state  space  Z  and  X  is  from  a  p  dimensional  parameter 
space.  In  general,  F  represents  some  boundary  value  problem  and  hence  the 
state  space  Z  is  infinite  dimensional.  Let  X  -  ZxA  be  the  combined  space, 
then  the  set  of  solutions  to  (1)  in  DcX, 

M  -  {(z,X)  £  D  |  F(z,X)  -  0)  (2) 

defines  a  surface  or  manifold  of  dimension  p  in  X.  M  is  called  the 
equilibrium  manifold. 

The  parameters  X  and  states  z  may  also  be  required  to  satisfy  a  differen¬ 
tial  equation  of  the  form 

A(x)x  -  G(x)  (3) 

where  x  ■  (z.X).  We  then  interpret  (3)  as  a  differential  equation  on  the  man¬ 
ifold  (2),  DEM  for  short.  Equations  (1)  and  (3)  together, 

f  F(x)  -  0 
U(x)x  -  G(x) 

form  what  is  called  a  differential -algebraic  equation  (DAE) . 

Two  applications  of  interest  to  the  investigators  in  which  DEM's  occur 

are : 

(i)  Incompressible  Fluid  Flow.  The  continuity  equation  is  an  algebraic  con¬ 
straint  of  the  form  (1)  and  the  time-dependent  Navier-Stokes  equations 
are  the  differential  equation. 

(ii)  Punch  Stretching  of  Sheet  Metal .  The  principle  of  virtual  work  provides 
a  force  equilibrium  equation  which  defines  an  equilibrium  manifold  upon 
which  one  seeks  solutions  to  differential  constitutive  laws  of  the  form 
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In  [20],  a  new  numerical  method  was  presented  for  computer  simulations  of 
punch  stretching  of  sheet  metal.  Most  current  approaches  to  finite  element 
modeling  of  large  deformation,  elastic-plastic  sheet  metal  forming  use  a  rate 
form  of  the  equilibrium  equations  and  then  must  correct  at  each  time  step  to 
insure  that  equilibrium  is  satisfied.  Such  methods  are  referred  to  as  incre¬ 
mental  methods .  The  new  method,  a  DEM  approach  discretizes  the  more  fundamen¬ 
tal  equilibrium  equations  in  non-rate  form  and  insures  equilibrium  of  forces 
at  each  time  step.  Formulating  the  problem  as  a  DEM  or  DAF  'Iso  allowed  for 
solution  of  the  discretized  system  using  off-the-shelf  software  such  as  LSODI . 
Numerical  experimentation  indicated  that  the  DEM  approach  was  computationally 
much  more  efficient  than  the  incremental  approach. 


Stability  of  Viscous  Incompressible  Flows 


The  problem  of  determining  sufficient  conditions  for  the  flow  of  a 
viscous  incompressible  fluid  to  be  stable  under  arbitrary  disturbances  was 
examined.  This  problem  is  of  importance  in  the  study  of  turbulence  and  the 
transition  which  occupies  a  region  of  space  and  subject  to  a  prescribed  velo¬ 
city  distribution  on  the  boundary,  will  alter  radically  or  only  slightly  in 
nature  if  it  is  perturbed  at  some  instance. 


The  question  of  stability  can  be  addressed  by  either  standard  linear  per¬ 
turbation  techniques  or  by  the  energy  method;  the  latter  is  chosen  in  this 
work.  Although  the  great  majority  of  stability  calculations  use  linear  sta¬ 
bility  analysis,  the  method  has  the  drawback  that  it  allows  only  perturbations 
which  are  infinitesimal  in  magnitude.  This  rules  out  perturbations  of  finite 
size  and  hence  cannot  give  accurate  information  in  many  cases.  The  energy 
method  allows  arbitrary  disturbances  but  its  shortcoming  is  that  the 
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disturbances  do  not  necessarily  satisfy  the  Navier- Stokes  equations,  and  thus 
the  stability  criterion  will  be  more  restrictive  than  in  the  actual  physical 
situation.  However,  the  energy  stability  analysis  is  based  on  the  Navier- 
Stokes  equations  and  is  nonlinear  in  nature  due  to  the  fact  that  no  lineariza¬ 
tions  of  the  equations  are  done.  The  method  is  mathematically  rigorous  and 
does  give  insight  into  the  physical  situation. 


The  question  of  energy  stability  of  a  flow  can  be  formulated  as  a  linear 
generalized  partial  differential  eigenvalue  problem  even  though  the  analysis 
is  based  on  the  nonlinear  Navier- Stokes  equations.  Essentially,  the  procedure 
is  to  obtain  an  equation  for  dK/dt  where  K  is  the  kinetic  energy  of  the  dis¬ 
turbance,  u,  and  then  to  determine  conditions  which  guarantee  that  the  kinetic 
energy  tends  to  zero  as  time  increases.  Using  standard  eigenvalue  problem 
where  stability  is  governed  by  the  dominant  eigenvalue.  Specifically,  we  have 
that  the  flow  is  stable  for  Reynolds  number  less  than  1/A  where  A  is  the  dom¬ 
inant  eigenvalue  of  the  problem 


Au  -  grad  4>  “  A  u  D 


div  u  -  0  , 


(1) 


u  -  0  on  the  boundary 

Here  D  is  the  deformation  tensor  of  the  unperturbed  flow. 


A  finite  element  method  is  used  to  approximate  the  dominant  eigenvalue  of 


(1).  In  particular,  the  weak  form  considered  is  to  find  nonzero  u  t  H  , 

—  — o 


4>  «  L  and  A  t  R  such  that 


J  Vu:Vv 


<f>  div  y  -  A  J  uDv 

*div  u  -  0 


for  all  v  c  H 
—  —  o 


for  all  x  f  L 


(2) 
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p 


A1. 

>  i 


2  ii 

where  L  is  the  space  of  all  functions  which  are  square  integrable  and  H  ,  H 

denote  the  usual  Sobolev  spaces.  To  approximate  the  solution  to  the  weak 

h  1  k  2 

problem  (2),  finite  dimensional  subspaces  V  c  Hx  and  w1  c  L  are  chosen  which 

—  — o 

depend  upon  a  parameter  0<h<l  tending  to  zero.  The  approximate  problem  is 
defined  analogous  to  (2)  where  the  solution  is  sought  in  the  finite  dimen¬ 
sional  subspaces.  Once  bases  for  and  W*1  are  chosen,  the  approximate  prob¬ 
lem  is  equivalent  to  an  algebraic  generalized  eigenvalue  problem.  An  estimate 
for  the  error  in  A  and  its  Galerkin  approximation  is  given  in  [A14] . 

As  proposed,  a  code  was  developed  which  uses  a  mixed  finite  element 
method  for  approximating  the  dominant  eigenvalue  of  (1).  The  program  was  used 
to  determine  a  range  of  Reynolds  numbers  for  which  the  flow  is  guaranteed  to 
be  stable  for  the  examples  of  plane  shear  flow  and  Poisseuille  flow. 

The  first  example  is  the  simple  case  of  flow  in  a  channel  of  width  0<y<d 
where  the  initial  velocity  is  given  by  the  vector  (ky,0),  the  deformation  ten¬ 
sor  is  given  by  -  0  i  -  j  and  -  .5k  for  i  *  j  ,  i,J  -  1,2  the  Reynolds 
2 

number  is  kd  /A.  The  channel  is  assumed  infinite  in  length.  The  computed 
Reynolds  number  for  a  channel  of  length  L  is  given  below. 

1/L  1  1/2  1/3  1/4  1/5 

Re  No.  289.87  201.42  186.975  182.11  180.01 

Using  the  above  calculations  the  extrapolated  value  at  1/L  -  0  is  177.4  which 
is  in  good  agreement  with  the  value  of  177  published  by  Orr. 

The  second  problem  of  determining  the  stability  of  Poiseuille  flow  in  a 
pipe  is  of  more  interest  physically  and  has  been  studied  by  many  authors.  For 
this  example  the  formulation  (1)  makes  sense  in  an  unbounded  domain  when  the 
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solutions  are  single-valued  in  6  and  periodic  in  z.  To  this  end,  the  solu¬ 
tion  is  assumed  of  the  form  u(r,*,z)  -  U(r,a,0)  ei(dZ  +  where  0  is  an 
integer.  This  form  is  substituted  into  (1)  and  the  following  system  results 


Lu  -  i  — ^  v  -  4>  ~  -Arw 

2  r 

r 


Lv  +  i(20/r  u  -  P/r4>)  -  0 


Lw  -  i a<f>  -  -Aru 


r  ji  <ru)  +  i(f  v  +  dw)  “  0 


u,v,w  bounded  at  r  -  0 


where  U  -  (u,v,w),  Lu  - 


u-v-w-Oatr-1 

-  (ru)  .  («!±1  .  A  ; 

* 


Lu  —  Lu  +  u/r 


Note  that  these  equations  form  a  complex  one -point  boundary  value  problem 
with  a  singularity  at  the  origin.  The  weak  formulation  must  incorporate  an 
appropriate  boundary  condition  for  the  velocity  at  r  -  0.  The  particular  weak 
from  considered  is  the  following:  Find  u.v.weH^  4>  e  H2>  A  e  R  such  that 


JUr^rr  dr  +  J  Clr  *dr  +  v£dr  "  J^~(rOdr  -  Ajw£r2dr  V  (  (  Hj 
|Vrrdr  +  J  Clr?dr  "  iJ^u^dr  +  ij^fdr  *  0 

JwrCrrdr  +  |  c2r^dr  +  iQJWrdr  -  *Jufr2dr  V  f  t  ^ 

-[^(ru)xdr  -  iJ/Sv^dr  -  ijaw*rdr  -  0  V  x  e  H2 


v  X  e  H, 


2  2  2 

where  C^-0  +l+or  •  c2  -  Cj  -  1 

. 

for  appropriate  spaces  and  H2 . 
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With  this  weak  formulation  the  condition  imposed  on  the  velocity  at  r  -  0 
is  ru^  -  rv  -  rw^  -  0,  which  is  a  natural  boundary  condition.  This  problem 
was  discretized  using  piecewise  linear  elements.  The  results  obtained  agree 
with  those  of  Joseph  and  Carmi  [A15J.  However,  their  numerical  calculations 
were  unnecessarily  complicated.  For  example,  different  techniques  for  the 
various  cases  such  as  a  -  0  and  a  *  0  had  to  be  employed  as  well  as  using  Fro- 
benius  series  as  starting  values  near  the  origin.  Specifically,  the  value  of 
81.5  was  obtained  as  a  sure  limit  of  stability  and  corresponded  to  the  case 
a  -  1,  (i  -  1.  This  agreed  with  Joseph  and  Carmi ‘s  result  and  confirmed  the 
fact  that  the  value  of  R  -  180,  which  was  previously  believed  to  be  a  sure 
limit  of  stability,  is  incorrect.  This  work  is  discussed  in  a  paper  in 
|  preparation  [23], 
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